#load all libraries and data file
library(survey)
library(ggplot2)
load(file="UMAS0035_PoP.RData")
stars.pval <- function (p.value){
  unclass(symnum(p.value, corr = FALSE, na = FALSE, cutpoints = c(0,0.001, 0.01, 0.05, 0.1, 1), symbols = c("***", "**","*", "'", " ")))}
p.legend <- "*** p<0.001 ** p<0.01 * p<0.05 ' p<0.1"

#Comparison of before and after, entire sample
lk.df <-data.frame(matrix(0,nrow=11,ncol=6))
colnames(lk.df) <- c("comp","ttest","lkf.mean","lkf.se","sqf.mean","sqf.se")

pn$sq.ge4.1a1b.2b2c <- NA
pn$sq.ge4.1a1b.2b2c[(pn$treatment==8|pn$treatment==9)] <- "2B+2C"
pn$sq.ge4.1a1b.2b2c[(pn$treatment==1|pn$treatment==2)] <- "1A+1B"
pn.design <- svydesign(ids=~0,data=pn,weights=~weight)
pn.svyby <- svyby(~flutest.eq2,~sq.ge4.1a1b.2b2c,design=pn.design,na.rm=TRUE,FUN=svymean)
pn.svyttest <- svyttest(flutest.eq2~sq.ge4.1a1b.2b2c,design=pn.design)
lk.df[1,1] <- "Entire Sample"
lk.df[1,2] <- pn.svyttest$p.value
lk.df[1,3:6] <- as.numeric(c(pn.svyby[2,2:3],pn.svyby[1,2:3]))

#Comparison of before and after for those who did and did not support the strike
#If you didn't support the strike (sq.ge4=0), does your law knowledge change depending on if you were asked about it before (8+9) or after (1+2)?
#did support strike
pn$sq.ge4.eq1.1a1b.2b2c <- NA
pn$sq.ge4.eq1.1a1b.2b2c[(pn$treatment==8|pn$treatment==9)&pn$sq.ge4==1] <- "2B+2C"
pn$sq.ge4.eq1.1a1b.2b2c[(pn$treatment==1|pn$treatment==2)&pn$sq.ge4==1] <- "1A+1B"
pn.design <- svydesign(ids=~0,data=pn,weights=~weight)
pn.svyby <- svyby(~flutest.eq2,~sq.ge4.eq1.1a1b.2b2c,design=pn.design,na.rm=TRUE,FUN=svymean)
pn.svyttest <- svyttest(flutest.eq2~sq.ge4.eq1.1a1b.2b2c,design=pn.design)
lk.df[2,1] <- "Prefer Strike"
lk.df[2,2] <- pn.svyttest$p.value
lk.df[2,3:6] <- as.numeric(c(pn.svyby[2,2:3],pn.svyby[1,2:3]))

#didn't support strike
pn$sq.ge4.eq0.1a1b.2b2c <- NA
pn$sq.ge4.eq0.1a1b.2b2c[(pn$treatment==8|pn$treatment==9)&pn$sq.ge4==0] <- "2B+2C"
pn$sq.ge4.eq0.1a1b.2b2c[(pn$treatment==1|pn$treatment==2)&pn$sq.ge4==0] <- "1A+1B"
pn.design <- svydesign(ids=~0,data=pn,weights=~weight)
pn.svyby <- svyby(~flutest.eq2,~sq.ge4.eq0.1a1b.2b2c,design=pn.design,na.rm=TRUE,FUN=svymean)
pn.svyttest <- svyttest(flutest.eq2~sq.ge4.eq0.1a1b.2b2c,design=pn.design)
lk.df[3,1] <- "Prefer Ground War"
lk.df[3,2] <- pn.svyttest$p.value
lk.df[3,3:6] <- as.numeric(c(pn.svyby[2,2:3],pn.svyby[1,2:3]))

#liberals
pn$ideo5.le2.1a1b.2b2c <- NA
pn$ideo5.le2.1a1b.2b2c[(pn$treatment==8|pn$treatment==9)&(pn$ideo5==1|pn$ideo5==2)] <- "2B+2C"
pn$ideo5.le2.1a1b.2b2c[(pn$treatment==1|pn$treatment==2)&(pn$ideo5==1|pn$ideo5==2)] <- "1A+1B"
pn.design <- svydesign(ids=~0,data=pn,weights=~weight)
pn.svyby <- svyby(~flutest.eq2,~ideo5.le2.1a1b.2b2c,design=pn.design,na.rm=TRUE,FUN=svymean)
pn.svyttest <- svyttest(flutest.eq2~ideo5.le2.1a1b.2b2c,design=pn.design)
lk.df[4,1] <- "Liberal"
lk.df[4,2] <- pn.svyttest$p.value
lk.df[4,3:6] <- as.numeric(c(pn.svyby[2,2:3],pn.svyby[1,2:3]))

#conservatives
pn$ideo5.ge4.1a1b.2b2c <- NA
pn$ideo5.ge4.1a1b.2b2c[(pn$treatment==8|pn$treatment==9)&(pn$ideo5==4|pn$ideo5==5)] <- "2B+2C"
pn$ideo5.ge4.1a1b.2b2c[(pn$treatment==1|pn$treatment==2)&(pn$ideo5==4|pn$ideo5==5)] <- "1A+1B"
pn.design <- svydesign(ids=~0,data=pn,weights=~weight)
pn.svyby <- svyby(~flutest.eq2,~ideo5.ge4.1a1b.2b2c,design=pn.design,na.rm=TRUE,FUN=svymean)
pn.svyttest <- svyttest(flutest.eq2~ideo5.ge4.1a1b.2b2c,design=pn.design)
lk.df[5,1] <- "Conservatives"
lk.df[5,2] <- pn.svyttest$p.value
lk.df[5,3:6] <- as.numeric(c(pn.svyby[2,2:3],pn.svyby[1,2:3]))

#registered to vote
pn$votereg.eq1.1a1b.2b2c <- NA
pn$votereg.eq1.1a1b.2b2c[(pn$treatment==8|pn$treatment==9)&(pn$votereg==1)] <- "2B+2C"
pn$votereg.eq1.1a1b.2b2c[(pn$treatment==1|pn$treatment==2)&(pn$votereg==1)] <- "1A+1B"
pn.design <- svydesign(ids=~0,data=pn,weights=~weight)
pn.svyby <- svyby(~flutest.eq2,~votereg.eq1.1a1b.2b2c,design=pn.design,na.rm=TRUE,FUN=svymean)
pn.svyttest <- svyttest(flutest.eq2~votereg.eq1.1a1b.2b2c,design=pn.design)
lk.df[6,1] <- "Registered to Vote"
lk.df[6,2] <- pn.svyttest$p.value
lk.df[6,3:6] <- as.numeric(c(pn.svyby[2,2:3],pn.svyby[1,2:3]))

#not registered to vote
pn$votereg.eq2.1a1b.2b2c <- NA
pn$votereg.eq2.1a1b.2b2c[(pn$treatment==8|pn$treatment==9)&(pn$votereg==2)] <- "2B+2C"
pn$votereg.eq2.1a1b.2b2c[(pn$treatment==1|pn$treatment==2)&(pn$votereg==2)] <- "1A+1B"
pn.design <- svydesign(ids=~0,data=pn,weights=~weight)
pn.svyby <- svyby(~flutest.eq2,~votereg.eq2.1a1b.2b2c,design=pn.design,na.rm=TRUE,FUN=svymean)
pn.svyttest <- svyttest(flutest.eq2~votereg.eq2.1a1b.2b2c,design=pn.design)
lk.df[7,1] <- "Not Registered to Vote"
lk.df[7,2] <- pn.svyttest$p.value
lk.df[7,3:6] <- as.numeric(c(pn.svyby[2,2:3],pn.svyby[1,2:3]))

#Just conservatives who opposed
pn$sq.ge4.eq0.ideo5.ge4.1a1b.2b2c <- NA
pn$sq.ge4.eq0.ideo5.ge4.1a1b.2b2c[(pn$treatment==8|pn$treatment==9)&(pn$ideo5==4|pn$ideo5==5)&pn$sq.ge4==0] <- "2B+2C"
pn$sq.ge4.eq0.ideo5.ge4.1a1b.2b2c[(pn$treatment==1|pn$treatment==2)&(pn$ideo5==4|pn$ideo5==5)&pn$sq.ge4==0] <- "1A+1B"
pn.design <- svydesign(ids=~0,data=pn,weights=~weight)
pn.svyby <- svyby(~flutest.eq2,~sq.ge4.eq0.ideo5.ge4.1a1b.2b2c,design=pn.design,na.rm=TRUE,FUN=svymean)
pn.svyttest <- svyttest(flutest.eq2~sq.ge4.eq0.ideo5.ge4.1a1b.2b2c,design=pn.design)
lk.df[8,1] <- "Conservative, Ground War"
lk.df[8,2] <- pn.svyttest$p.value
lk.df[8,3:6] <- as.numeric(c(pn.svyby[2,2:3],pn.svyby[1,2:3]))

#Just conservatives who supported
pn$sq.ge4.eq1.ideo5.ge4.1a1b.2b2c <- NA
pn$sq.ge4.eq1.ideo5.ge4.1a1b.2b2c[(pn$treatment==8|pn$treatment==9)&(pn$ideo5==4|pn$ideo5==5)&pn$sq.ge4==1] <- "2B+2C"
pn$sq.ge4.eq1.ideo5.ge4.1a1b.2b2c[(pn$treatment==1|pn$treatment==2)&(pn$ideo5==4|pn$ideo5==5)&pn$sq.ge4==1] <- "1A+1B"
pn.design <- svydesign(ids=~0,data=pn,weights=~weight)
pn.svyby <- svyby(~flutest.eq2,~sq.ge4.eq1.ideo5.ge4.1a1b.2b2c,design=pn.design,na.rm=TRUE,FUN=svymean)
pn.svyttest <- svyttest(flutest.eq2~sq.ge4.eq1.ideo5.ge4.1a1b.2b2c,design=pn.design)
lk.df[9,1] <- "Conservative, Strike"
lk.df[9,2] <- pn.svyttest$p.value
lk.df[9,3:6] <- as.numeric(c(pn.svyby[2,2:3],pn.svyby[1,2:3]))

#Just liberals who opposed
pn$sq.ge4.eq0.ideo5.le2.1a1b.2b2c <- NA
pn$sq.ge4.eq0.ideo5.ge2.1a1b.2b2c[(pn$treatment==8|pn$treatment==9)&(pn$ideo5==2|pn$ideo5==1)&pn$sq.ge4==0] <- "2B+2C"
pn$sq.ge4.eq0.ideo5.ge2.1a1b.2b2c[(pn$treatment==1|pn$treatment==2)&(pn$ideo5==2|pn$ideo5==1)&pn$sq.ge4==0] <- "1A+1B"
pn.design <- svydesign(ids=~0,data=pn,weights=~weight)
pn.svyby <- svyby(~flutest.eq2,~sq.ge4.eq0.ideo5.ge2.1a1b.2b2c,design=pn.design,na.rm=TRUE,FUN=svymean)
pn.svyttest <- svyttest(flutest.eq2~sq.ge4.eq0.ideo5.ge2.1a1b.2b2c,design=pn.design)
lk.df[10,1] <- "Liberal, Ground War"
lk.df[10,2] <- pn.svyttest$p.value
lk.df[10,3:6] <- as.numeric(c(pn.svyby[2,2:3],pn.svyby[1,2:3]))

#Just liberals who supported
pn$sq.ge4.eq1.ideo5.le2.1a1b.2b2c <- NA
pn$sq.ge4.eq1.ideo5.le2.1a1b.2b2c[(pn$treatment==8|pn$treatment==9)&(pn$ideo5==2|pn$ideo5==1)&pn$sq.ge4==1] <- "2B+2C"
pn$sq.ge4.eq1.ideo5.le2.1a1b.2b2c[(pn$treatment==1|pn$treatment==2)&(pn$ideo5==2|pn$ideo5==1)&pn$sq.ge4==1] <- "1A+1B"
pn.design <- svydesign(ids=~0,data=pn,weights=~weight)
pn.svyby <- svyby(~flutest.eq2,~sq.ge4.eq1.ideo5.le2.1a1b.2b2c,design=pn.design,na.rm=TRUE,FUN=svymean)
pn.svyttest <- svyttest(flutest.eq2~sq.ge4.eq1.ideo5.le2.1a1b.2b2c,design=pn.design)
lk.df[11,1] <- "Liberal, Strike"
lk.df[11,2] <- pn.svyttest$p.value
lk.df[11,3:6] <- as.numeric(c(pn.svyby[2,2:3],pn.svyby[1,2:3]))

write.csv(lk.df,file="lk.df.csv")

#plots
lk2.df <- lk.df[1:11,c(1,2,5,6)]
dimnames(lk2.df)[[2]][3:4] <- dimnames(lk.df)[[2]][3:4]

lkplot.df <- data.frame(rbind(lk.df[1:11,1:4],lk2.df))
dimnames(lkplot.df)[[2]][1:4] <- c("x","ttest","dv","se")
lkplot.df$ord <- paste0(lkplot.df$x,c(rep("1",11),rep("2",11)))
lkplot.df$comp <- c(rep("Before Iran Scenario Question",11),rep("After Iran Scenario Question",11))
lkplot.df$p.stars <- paste0("P=",signif(lkplot.df$ttest,2),stars.pval(lkplot.df$ttest))
lkplot.df$x.stars <- paste0(lkplot.df$x,stars.pval(lkplot.df$ttest))
lkplot.df$x <- factor(lkplot.df$x,levels=lkplot.df$x[1:11])
lkplot.df$x.stars <- factor(lkplot.df$x.stars,levels=lkplot.df$x.stars[1:11])
lkplot.df$comp <- factor(lkplot.df$comp,levels=c("Before Iran Scenario Question","After Iran Scenario Question"))

#figure 1
lkplot.df <- lkplot.df[c(1:3,12:14),]
ggplot(data.frame(lkplot.df),aes(x=x.stars,y=dv,fill=comp))+
  geom_col(position="dodge")+
  geom_errorbar(position=position_dodge(0.9),aes(ymin=dv-se,ymax=dv+se),width=0.4)+
  #  geom_text(aes(x=x,y=1.05,label=p.stars))+
  geom_text(position=position_dodge(0.9),aes(y=dv+0.01,label=paste0(signif(dv*100,2),"%")))+ #plot percentages
  labs(x=p.legend,y="Law Knowledge Correct",fill="If Asked...")+
  scale_y_continuous(breaks=c(0,.25,.5,.75,1),labels = scales::percent_format(accuracy=1),limits=c(0,1))+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
#  guides(fill=FALSE) #no legend
filename <- "Figure1.t1a1b.2b2c-3.pdf"
dev.copy(pdf,filename,width=6,height=5)
dev.off()

#verify that treatment*strike preference is insignificant: footnote 11
pn$t1a1b.2b2c <- NA
pn$t1a1b.2b2c[(pn$treatment==8|pn$treatment==9)] <- 0
pn$t1a1b.2b2c[(pn$treatment==1|pn$treatment==2)] <- 1
pn.design <- svydesign(ids=~0,data=pn,weights=~weight)
summary(svyglm(flutest.eq2 ~ t1a1b.2b2c + sq.ge4 + t1a1b.2b2c*sq.ge4, design = pn.design, family=quasibinomial()))
#and other demographic variables
summary(svyglm(flutest.eq2 ~ t1a1b.2b2c + gender.bin + t1a1b.2b2c*gender.bin, design = pn.design, family=quasibinomial()))
summary(svyglm(flutest.eq2 ~ t1a1b.2b2c + ideo5 + t1a1b.2b2c*ideo5, design = pn.design, family=quasibinomial()))
summary(svyglm(flutest.eq2 ~ t1a1b.2b2c + pid7 + t1a1b.2b2c*pid7, design = pn.design, family=quasibinomial()))
summary(svyglm(flutest.eq2 ~ t1a1b.2b2c + presvote16post.bin + t1a1b.2b2c*presvote16post.bin, design = pn.design, family=quasibinomial()))
summary(svyglm(flutest.eq2 ~ t1a1b.2b2c + educ-1 + t1a1b.2b2c*(educ-1), design = pn.design, family=quasibinomial()))

#Figure 2
pn$t1a1b.1c <- NA
pn$t1a1b.1c[pn$treat.nam=="1A"|pn$treat.nam=="1B"] <- 1
pn$t1a1b.1c[pn$treat.nam=="1C"] <- 2
pn.design <- svydesign(ids=~0,data=pn,weights=~weight) #set up survey design
pn.svyttest <- svyttest(sq.gt5lt2~t1a1b.1c,pn.design,na.rm=TRUE) #sig!
pn.svyby <- svyby(~sq.gt5lt2,~t1a1b.1c,design=pn.design,na.rm=TRUE,FUN=svymean)
pn.svyby$p.value <- pn.svyttest$p.value
names(pn.svyby)[names(pn.svyby)=="sq.gt5lt2"] <- "dv"
filename <- "Figure2.gt5lt2.t1a1b.tc.pdf"
dv.nam <- "Strike Support (strong preferences only)"
pn.svyby$p.stars <- paste0("P=",signif(pn.svyby$p.value,2),stars.pval(pn.svyby$p.value))
pn.svyby$label <- factor(c("No Survey Information","Survey Information"))
ggplot(pn.svyby,aes(x=t1a1b.1c,y=dv,fill=label))+
  geom_col()+
  #  ggtitle("Treatments 1-2 Strike Support v Treatment 7 Strike Support (after sentiment)")+
  geom_text(aes(x=1.5,y=max(dv+se)+0.05,label=p.stars))+
  #scale_fill_manual(values=rainbow(4)[c(4,3)])+
  geom_errorbar(aes(ymin=dv-se,ymax=dv+se),width=0.4)+
  ylab(dv.nam)+
  geom_text(aes(y=dv,label=paste0(signif(dv*100,2),"%")))+ #plot percentages
  scale_y_continuous(breaks=c(0,.25,.5,.75,1),labels = scales::percent_format(accuracy=1),limits=c(0,1))+
  scale_x_continuous(breaks=c(1,2),name=element_blank(),labels=pn.svyby$label)+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  guides(fill=FALSE) #no legend
dev.copy(pdf,filename,width=3.5,height=5)
dev.off()

